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MECHANICAL BALANCE LAWS FOR FULLY NONLINEAR 
AND WEAKLY DISPERSIVE WATER WAVES 

HENRIK KALISCH, ZAHRA KHORSAND, AND DIMITRIOS MITSOTAKIS 


Abstract. The Serre-Green-Naghdi system is a coupled, fully nonlinear sys¬ 
tem of dispersive evolution equations which approximates the full water wave 
problem. The system is known to describe accurately the wave motion at the 
surface of an incompressible inviscid fluid in the case when the fluid flow is ir- 
rotational and two-dimensional. The system is an extension of the well known 
shallow-water system to the situation where the waves are long, but not so 
long that dispersive effects can be neglected. 

In the current work, the focus is on deriving mass, momentum and energy 
densities and fluxes associated with the Serre-Green-Naghdi system. These 
quantities arise from imposing balance equations of the same asymptotic or¬ 
der as the evolution equations. In the case of an even bed, the conservation 
equations are satisfied exactly by the solutions of the Serre-Green-Naghdi sys¬ 
tem. The case of variable bathymetry is more complicated, with mass and mo¬ 
mentum conservation satisfied exactly, and energy conservation satisfied only 
in a global sense. In all cases, the quantities found here reduce correctly to 
the corresponding counterparts in both the Boussinesq and the shallow-water 
scaling. 

One consequence of the present analysis is that the energy loss appear¬ 
ing in the shallow-water theory of undular bores is fully compensated by the 
emergence of oscillations behind the bore front. The situation is analyzed 
numerically by approximating solutions of the Serre-Green-Naghdi equations 
using a finite-element discretization coupled with an adaptive Runge-Kutta 
time integration scheme, and it is found that the energy is indeed conserved 
nearly to machine precision. As a second application, the shoaling of solitary 
waves on a plane beach is analyzed. It appears that the Serre-Green-Naghdi 
equations are capable of predicting both the shape of the free surface and the 
evolution of kinetic and potential energy with good accuracy in the early stages 
of shoaling. 


1. Introduction 

In this paper we study mechanical balance laws for fully nonlinear and disper¬ 
sive shallow-water waves. In particular, the Serre-Green-Naghdi (SGN) system of 
equations with variable bathymetry is considered. This system was originally de¬ 
rived for one-dimensional waves over a horizontal bottom in 1953 by F. Serre Hi- 
Several years later, the same system was rederived by Su and Gardner [3b In 1976, 
Green and Naghdi [4] derived a two-dimensional fully nonlinear and weakly dis¬ 
persive system for an uneven bottom while in one spatial dimension Seabra-Santos 
et al. [5] derived the generalization of the Serre system with variable bathymetry. 
Lannes and Bonneton derived several other systems including the SGN equations 
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using a new formulation of the water wave problem, [6j. For more information and 
generalizations of the SGN equations we refer to [7] and the references therein, 
while we refer to the paper by Barthelemy [B] for an extensive review. 

The Serre-Green-Naghdi (SGN) system and several variants of it are extensively 
used in coastal modeling pm eh [ 3 . In the present contribution, the focus is on 
one aspect of the equations which has not received much attention so far, namely 
the derivation and use of associated mechanical balance equations, and in partic¬ 
ular a differential energy balance equation. While it is known that the equations 
admit four local conservation equations if the bed is even [12], it appears that the 
connection to mechanical balance laws of the original Euler equations has not been 
firmly established so far. Here we show that the first three conservation laws of the 
Serre-Green-Naghdi (SGN) equations arise as approximations of mechanical bal¬ 
ance laws in the context of the Euler equations, both in the case of even beds, and 
in the case of nontrivial bathymetry. The fourth conservation law has been shown 
to arise from a kinematic identity similar to Kelvin’s circulation theorem 1131 - 

Let us first review some modeling issues regarding the Serre-Green-Naghdi (SGN) 
system. Suppose a denotes a typical amplitude, and l a typical wavelength of a 
wavefield under study. Suppose also that bo represents the average water depth. In 
order to be a valid description of such a situation, the SGN equations require the 
shallow water condition, /3 = b^/l 2 -C 1. In contrast, the range of validity of the 
weakly nonlinear and weakly dispersive Boussinesq equations is limited to waves 
with small amplitude and large wavelength, i.e. a = a/bo 1 and /3 <C 1. In 
this scaling regime, one also finds the weakly nonlinear, fully dispersive Whitham 
equation (annus]. 

The SGN equations can be derived by depth-averaging the Euler equations and 
truncating the resulting set of equations at 0(/3 2 ) without making any assumptions 
on the order of a , other than a < 0(1). 

In their dimensionless and scaled form the SGN equations can be written as 

Vt + [hu] x = 0 , (la) 

Ut + uu x + gr] x + — [h 2 {^P + |S)] I + b x (^V + Q) = 0 , (lb) 

with V = h [u 2 — u x t — uu xx ] and Q = b x (ut + uu x ) + b xx u 2 , x € R, t > 0, along 
with the initial conditions h(x,0) = ho(x), u{x, 0) = uo{x). Here, ij = r)(x,t) is the 
free surface displacement, while 

h = r]-b , (2) 

denotes the total fluid depth. The unknown u = u(x,t) is the depth-averaged 
horizontal velocity, and rjo, uq are given real functions, such that rjo — b > 0 for all 
x £ R. In these variables, the location of the horizontal bottom is given by z = b 
(cf. Figure 1). For a review of the derivation and the basic properties of this system 
we also refer to mm- 

In the case of small-amplitude waves, i.e. if (3 ~ a, the SGN equations reduce to 
Peregrine’s system HU- On the other hand, in the case of very long waves, i.e. /3 —>• 
0, the dispersive terms disappear, and the system reduces to the nondispersive 
shallow water equations. 

The SGN system for waves over a flat bottom possesses solitary and cnoidal wave 
solutions given in closed form. For example, the solitary wave with speed c s can be 
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written as 


MO = h s {x,t) = o 0 + aisech 2 (/v s £), 


u s {0 = U s (x,t) = Ca ( 1 - 


ao 

MO 


(3a) 

(3b) 


where £ = a; — c s t, K s = y/3a\/ Aa^c 2 , c s = y/ao~ FaT, and ao > 0 and ai > 0. 
For more information about the solitary and cnoidal waves and their dynamical 
properties we refer to emails uses Eg. 

It is important to note that the SGN system has a Hamiltonian structure, even 
in the case of two-dimensional waves over an uneven bed cf. US 121 H2123- Specif¬ 
ically, any solution ( h,u ) of (jTJ) conserves the Hamiltonian functional 

1 r°° 

'H{t) = - / gr 1 2 + hu 2 + h [h x b x + \hb xx + b 2 x ] u 2 - | [h 3 u x \ x u dx , (4) 

in the sense that dH(t)/dt = 0. Note however that (Hal) . (Ilbl) are recovered only if a 
non-canonical symplectic structure matrix is used. While in many simplified models 
equations, the Hamiltonian functions does not represent the mechanical energy of 
the wave problem [25], in the case of SGN, the Hamiltonian does represent the 
approximate total energy of the wave system. Thus the Hamiltonian can be written 
in the form 


/ OO 

E(x,t) 

-OO 


dx 


where the integrand 
1 


E = - (grf + hu 2 + h [h x b x 


2 tib xx 


u 2 


[h 3 Ux ] x u) 


= 0 


(5) 


is the depth-integrated energy density. In the present paper, we also identify a local 
depth-integrated energy flux qe, such that an equation of the form 

dE dqs 
dt dx 

is satisfied approximately. The procedure of finding the quantities E and qE follows 
a similar outline as the derivations in : 25j for a class of Boussinesq systems and m 
for the KdV equation. Expressions for energy functionals associated to Boussinesq 
systems have also been developed in [28 j. 

The analytical results are put to use in the study of undular bores. It is well 
known that the shallow-water theory for bores predicts an energy loss [25]. In 
an undular bore, the energy is thought to be disseminated through an increasing 
number of oscillations behind the bore, and the traditional point of view is that 
dissipation must also have an effect here [30] [31]. However, recent studies [32] have 
shown that if dispersion is included into the model equations, then the energy loss 
experienced by an undular bore can be accounted for without making appeal to 
dissipative mechanisms. 

Indeed, it was argued in ,32, 33 ] that the energy loss in an undular bore could 
be explained wholly within the realm of conservative dynamics by investigating 
a higher-order dispersive system, and monitoring the associated energy functional. 
However, there was a technical problem in the analysis in these works, as the energy 
functional was not the same as the one required by the more in-depth analysis in 
[25] , On the other hand, the energy functional found in [2B[ did not reduce to 
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the shallow-water theory in the correct way. In the current contribution, it is 
our purpose to remedy this situation by using the SGN system which reduces in 
the correct way to the shallow-water equations, and also features exact energy 
conservation in the case of a flat bed. 

The numerical method to be used is a standard Galerkin / Finite Element 
Method (FEM) for the SGN equations with reflective boundary conditions extend¬ 
ing the numerical method presented in |34j . For the sake of completeness we men¬ 
tion that there are several numerical methods applied to boundary value problems 
of the SGN equations. For example finite volume 12511351133, finite differences 
PH [38] [39] > spectral [36] and Galerkin methods [34] [40] . 

The paper is organized as follows: A review of the derivation of the SGN equa¬ 
tions based on um is presented in Section [2] The derivation of the mass, mo¬ 
mentum and energy balance laws in the asymptotic order of the SGN equations 
is presented in Section [3] Applications to undular bores and solitary waves are 
discussed in Section 0] The numerical method to be used used in this paper is 
presented briefly in 1X1 

2. The SGN equations over a variable bed 

Before introducing the balance laws for the SGN equations, we briefly review the 
derivation of the SGN equations from the Euler equations following the work [8], 
but in the case of a general bathymetry. This well known derivation is included here 
to set the stage for the development of the approximate mechanical balance laws in 
the next section. We consider an inviscid and incompressible fluid, and assume that 
the fluid flow is irrotational and two-dimensional. Let ao be a typical amplitude, 
l a typical wavelength and bo a typical water depth. We perform the change of 
variables x = x/l, z = z/bo, t = Cot/l, which yields non-dimensional independent 
variables identified by tildes, where x represents the horizontal and z the vertical 
coordinate. The limiting long-wave speed is defined by cq = y/gb o, and g denotes the 
acceleration due to gravity. The non-dimensional velocity components are defined 
by u = u/aco , v = v/y/]3aco, where a = ao/bo and /? = b^/l 2 . Finally, the free 
surface deflection, bottom topography and pressure are non-dimensionalized by 
taking rj = rj/ao , b = b/bo, and p = p/ pgbo- 


Z 



Figure 1. The geometry of the problem 
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In non-dimensional variables, the free-surface problem is written as follows [4Tj : 


The momentum equations are 

aui + a 2 (u 2 )x + a 2 (uv)z = -pi , (6a) 

+ a 2 fiuvi + a 2 pvvs = —p; — 1 . (6b) 

The equation of continuity and the irrotationality are expressed by 

us + Vz = 0 , (7a) 

uz - Pvt = 0 . (7b) 


The boundary conditions at the free surface and at the bottom are given by 

v = pi + aufjx j at 5 = aij(x) , (8a) 

p = 0, at z = afj(x) , (8b) 

v = bxU, at z = b(x) . (8c) 


The first equation in the system m is obtained by integrating the equation of 
continuity over the total depth. The result is written in terms of the depth-averaged 
horizontal velocity 

1 f af > 

u = — u dz , (9) 

hJi 

in the form 


Vt + [hu]x = 0 . (10) 

Using the boundary conditions (15al) -(15cl). the continuity equation m and the 
depth-averaged momentum equation (lirall yields 


d f v ~ I" v 

ahui + a 2 huux + q 2 tt / ( u 2 — (w) 2 ) dz = — / pi dz . 

ox Jj, Ji 

Applying the Leibniz rule to the right-hand side of equation m yields 


( 11 ) 


d 

Pi dz = — (hp) - amp\ i=afj + bip\~ =l 


d 

= g=( h P) + h *P\~z=b ■ 

(12) 

The momentum equation (16bl) is rewritten as 


a/3T(x,z, t) = -l~Pz , 

(13) 

where 


r(x, z, t) = vi + auvi + avvz ■ 

(14) 

Integrating equation (fl3l) from z to ap yields 


f af) 

p{x,z,t) = (otfj - z)+ a/3 r(x,(,t)d(, 

J Z 

(15) 

and taking the mean value gives 


^ rCX-fj ncx.fi 

hp =-h 2 + a(3 / / T(x,£,t)d£ dz . 

2 Ji Jz 

(16) 
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Therefore, equation m can be written as 

pd[ af, ~ 

ux + auux + Vx + — / (z — b)T(x,z,t)dz 

h ox J~ b 

o ncxfj _ c\ potfj 

+ jbx / T(x,z,t)dz = ~ r — (u 2 - (u) 2 ) dz . 

h Jb h Ox Ji 

The non-dimensional velocity componens are given (cf. |16] 1 to first order by 


u(x, 5, t) = u(x, t) + 0(j3) , 

(17) 

~ ~ i — f)h 

v(x,z,t) = -(z- b{x))— + 0(13) . 

(18) 


As it was shown in m, we can expand the velocity components using Taylor 
series in the vertical coordinate around the bottom. Denoting by u b and v b , respec¬ 
tively, the horizontal and vertical velocities at the bottom, the bottom kinematic 
condition (18cl) imposes that v b = b^u b . In order to determine which terms should 
be kept to obtain an approximation for the velocity field, the incompressibility con¬ 
dition (I7al) must hold to the same order in (3 as the evolution equations. If the 
non-dimensional velocity components are given by 

u(x, z,t) = u b (x,t)+/3(z-b) (bxul + (pxU b )^j -^(z-b) 2 ul s +0(/3 2 ) , (19) 


v{x, z , t ) = bxU b + (z~b) + P(bx(u b bx)x + &!?>!)) 

- f (z - bf (bsculs. + (b- x u b - + (ibxu b )x)x ) + ^(5 - bfu + 0((3 2 ) , 

then the incompressibility condition (f7a|) holds to 0(p 2 ). Depth averaging (fT9|) 
gives 

u b = u- (biUx + (hu)ij + ^h 2 Uxx + 0((3 2 , a/3 2 ) . 

Thus the horizontal velocity is 


i(x, z,t) = u — (3 (bxiix + (bxu)z) ( ^ - (z - b) 


+ A } ^\{z-b) 2 \^i + 0{p 2 ,ap 2 ). ( 20 ) 


Taking squares in equation o 


u 2 (x, z,t) = u 2 — f3 (bxUxii + ( bxu)xuj (h — 2 (z — b)j + 

13 - (z - b) 2 ^j uUxx + 0(/3 2 ,a/3 2 ) . (21) 

Integrating equation m from b to afj and after some simplifications it follows that 


ncx.fi 

/ (u 2 - (u) 2 ) dz = 0(/3 2 ,a/3 2 ) 

Jb 


(22) 
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and that 


T(i, 3, t) = (3 - b) [cm? - u^ - auu xx \ + 

+ bx{uf + auu x ) + abz.xu 2 + 0{/3, a/3) 

Evaluating the integrals T dz and J^iz — b)T dz yields 

pafj i 

/ Tdz= -hV + hQ , 


and 


where 


and 


r ar > _ 1 _ 1 _ 

(.3 — b)T dz = -h 2 V + -h 2 Q 

k o 2 


V = h \au\ — u%i — auu xx \ 


(23) 

(24) 

(25) 

(26) 
(27) 


Q = b x (u t ~ + auu x ) + b xx u 2 . 

Finally we find the second equation of the system as 

ui + auu x + fjx + ■=■ -gZ | \ h 2 1 + f3b x ^5 “P + = 0{aj3 2 ) . 

By setting the right-hand side equal to zero, and writing the variables in dimen¬ 
sional form the system reads 


Vt + [hu] x = 0 , 

Ut + uu x + gr] x + — [h 2 (jfP + 5Q )] x + b x {\£P + Q) = 0 , 


(28a) 

(28b) 


where V = h [u 2 — u x t — uu xx \ and Q = b x (ut + uu x ) + b xx u 2 . 

In order to determine which terms should be kept for the velocity field at a 
certain order of approximation, the incompressibility condition (17al) can be used. 
Then, the dimensional form of the water particle velocities at any location (x, z) in 
the vertical plane become 


u = u 


v = — zu r 


h 2 z 2 


(29a) 

(29b) 


As it was mentioned before, system (Hal) and m reduces to the shallow water 
system when /3 —> 0 and to the classical Boussinesq system when j3 ~ a. 

An asymptotic expression for the pressure p(x,z,t) can be obtained by substi¬ 
tuting formula (1231) into m Such a formula was derived in [12] in the form 


_. „ „ _ a /3 r — — — 2 

pyx, z,t) = arj — z + - 1 —— [-u x i - auu xx + au x 


](' 


h 2 — (3 — 5 ) 2 ) 

+ a/3 (ab xx u 2 + ab x uu x + b x uA (afj — z ) + 0{a/3 2 ) . (30) 
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3. Mechanical balance laws for the SGN equations 

In this section we derive the mechanical balance laws such as the mass, momen¬ 
tum and energy conservation for the SGN equations extending the results related 
to some Boussinesq systems found in |26 . The balance laws consist of terms of the 
same asymptotic order as in the SGN equations. We start with the conservation of 
mass. 


3.1. Mass balance. We investigate the mass conservation properties of equations 
SaT) and (I28bl) . Our starting point is the total mass of the fluid contained in a 


control volume of unit width, bounded by the lateral sides of the interval [xi,X 2 ], 
and by the free surface and the bottom. This mass is given by 


M = 


pdzdx. 


(31) 


According to the principle of mass conservation and the fact that there is no mass 
flux through the bottom or the free surface, mass conservation can be considered 
in terms of the flow variables as follows: 


d 

dt 


pdzdx = 


pu(x, z, t) dz 


(32) 


In non-dimensional form this equation becomes 


d_ 

dt 


/•ar] 


dzdx = a 


u(x, z, t ) dz 


(33) 


Substituting the expression (BUI) for u and integrating with respect to z yields 
d 


— / hdx = a 
dt J 


uh 


+ 0(ap 2 ) 


(34) 


where h = afj — b denotes the nondimensional total depth. Dividing (1341) by X 2 — X\ 
and taking x? — X\ -A 0 then the mass balance equation is written as 

hf + ( auh)x = 0(a/3 2 ) . (35) 

Denoting the non-dimensional mass density by M = h and the non-dimensional 
mass flux by qm = auh, then the mass balance is 


dt ox 


(36) 


Using the scaling M = pb$M and qM = pboCoqM the dimensional forms of mass 
density and mass flux are M = ph and qM = puh respectively. Then the dimen¬ 
sional form of the mass balance is obtained by discarding the right-hand side of the 
scaled mass balance equation and using the unsealed quantities: 


dM dq M 
dt dx 


= 0 


(37) 


It is noted that the mass balance is satisfied exactly by the solutions of the SGN 
system. 

The expressions for mass density and the mass flux do not depend on the shape 
of the bottom topography, and in particular, they have the same form for both even 
and uneven beds. The dimensional form of (1351) coincides with analogous formulas of 
the shallow-water wave system and the classical Boussinesq system [2U]. While this 













MECHANICAL BALANCE LAWS FOR FULLY NONLINEAR AND WEAKLY DISPERSIVE WATER WAVES 


may be expected, it should be pointed out that in the case of other asymptotically 
equivalent systems, mass conservation may be satisfied only to the same order as 
the order of the equations, [25] , 

3.2. Momentum Balance. The total horizontal momentum of a fluid of constant 
density p contained in a control volume of the same type as in the previous section 
is 

fX2 fV 

pu dz dx . (38) 


X = 


rx 2 l-r, 

J x i J b 


Conservation of momentum implies that the rate of change of I is equal to the net 
influx of momentum through the boundaries plus the net force at the boundary of 
the control volume. Therefore, the conservation of momentum is written 


d 

dt 


pu dz dx = — / pb x dx + 


rV rV ' 

/ pu 2 (x, z) dz + / pdz 

Jb Jb 


Non-dimensionalization of this expression leads to 

pX2 rotf) rX2 


d 

a—= 
dt , 


udzdx = — 


Pbbx dx + 


.af) 

.af) 

a 2 

u 2 dz+ pdz 

Jb 

Jb 


where Pb denotes the pressure at the bottom Pb = h+ ^r[— u^— otuuxx + cm|]/i 2 + 
a/3{abxxU 2 + abxUUx + bxiiph. Substituting the values of u and p from equations 
m and ([T5D and integrating with respect to z yields 


d f X2 = ~ 

a— / uhdx = — 

dt J Xl 


P b bx dx 


a 2 u 2 h + y - + auiixx - a(iix) 2 ) 


y- h 2 (abxxU 2 + bx (auux + up) 


0(a/3 2 


Applying similar techniques used for the derivation of the mass balance equation 
we obtain the momentum balance equation in the form 


(®uh)_ + | a 2 u 2 h + y - -y^ 3 - a{ux) 2 ) J + 




[ abnu + bx{auux + y 

If the non-dimensional momentum density is defined by 
I = auh 

and the momentum flux plus pressure force is defined by 

- 2~2L d Otfd ~ o — —— — 2 , 

qi = a u h T —-y h (u 2t - + auuxx - ct{Ux) )+ 


= ~P b h + 0(a/3 2 ) 


(39) 

(40) 


~ ~ _ ~ _ _ 

—h 2 (abxxU 2 + bx(auux + up) 
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then the momentum balance equation can be written as 

f + § = -^ + °(“«- < 41 > 

Using the scaling / = pcobol and qj = pc^b^qj, the dimensional forms of the 
momentum density and momentum flux per unit span are given by 

/ = puh , (42) 


and 


qi = pu 2 h + !yh 2 - ^(u xt + uu xx - u 2 x )h 3 + 

+ b xx u 2 + b x (uu x + u t ))h 2 , (43) 

respectively. It turn out that the momentum conservation law is also satisfied 
exactly in the context of the SGN system. Indeed, if the momentum density is 
defined by (l42l) . the momentum flux plus pressure force is defined by (l43l) . and the 
pressure is defined by fUD, then solutions of the SGN system also satisfy exactly 
the equation 


d]_ dqi_ 

dt dx 


—b x P ■ 


(44) 


Note that if the bottom b = —bo is horizontal, then the last equation is homogeneous 
and does not depend on the pressure p. 

Taking /3 —> 0 in the momentum balance equation (1391) . and using dimensional 
variables and horizontal bottom b = —bo , the momentum density is unchanged, but 
the flux reduces to 


dT 


-2 l l P9 72 
= pu h + —h. 


(45) 


Thus it is plain that both the momentum density I and flux qj reduce correctly to 
the nonlinear shallow water approximation. In the case /3 ~ a and a flat bottom, the 
quantities for the momentum balance law are I = pu(bo+rf) and qi = pb^u 2 +^h 2 — 
f&o'iixtj which agree with the corresponding quantities of the classical Boussinesq 
system. 


3.3. Energy Balance. The total mechanical energy inside a control volume can 
be written as the sum of the kinetic and potential energy as 


£ = 


j-x 2 rri 

/ / {f (u 2 + v 2 ) + pgz) dzdx . 

J x i J b 


(46) 


The conservation energy can be expressed as 
(■x 2 r v 


d_ 

dt 


px 2 rr] 

/ / {%(u 2 + v 2 ) + pgz} dzdx = 

J x i J b 

" n 

/ {(%{u 2 + v 2 ) + pgz) u + uP} dz 

Jb 


( 47 ) 
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and in non-dimensional variables as 
d 
dt 


f-x 2 rotr) 

' x i Jb 


j^-(u 2 + (3v 2 ) + zj dzdx = 

r af > r 2 i 

J |^-(u 3 + (3v 2 u) + zu + puj dz 


• (48) 


By substituting the expressions G 3 , m and (l30l) for u, v and p respectively, the 
energy balance equation takes the form 


d 

dt 


[u 2 + (3b 2 u 2 ^ h - ^y-b £ h 2 uus-+ 

^-u 3 (l + (3b 2 ^ h + 7}h 2 u+ 


a 2 (3 73~2 h 2 ~~ \ 

— —h 3 ui H-b bh \ dx = 

6 2 / 


0/3 Pi ~ ~2l2 i ^ ft ~.~,2 j L i 


abuh - —bxUxU~h~ -\——uri 


a 9 a (3 ~ o — / — —— —o 

— uh~ — h uyu~i + auu^x ~ cm 


) 


a 2 (3 


[abxx'U 2 + abx{uux + h 2 


0(a(3 2 ). (49) 


The differential form of the energy balance equation is given by 

{ a ~ ( ~2 , oZ2 ~2\ I 0/2 ft U u2~~ , P u3~2 i^ 2 

+ fib s u Jh - —bxh uus + —/Cm 5 + — + 


bh 


+ ( —vf'h 4 - -(—b^.u' 5 + auh 2 + abuh — ^ l3 ~ / 


: 3 7 a 3 (3 ~2 -3 

2 dh + ^ 
a 3 (3 - 


«Pl3d 

— h »( 


u Z i + auuxx 


--auz. - 


OxU £ U 




abxxu 2 + bx(oiuux + Uf )) = 0(a/3 2 ) . (50) 


Considering the appropriate terms in the energy density and flux in (I48|l which 
are of order zero or one in the differential energy balance (l49l) . we find that the 
non-dimensional energy density is 


E = ^-(u 2 + (3blu 2 )h~ bxh 2 uux + + bh , (51) 

while the non-dimensional energy flux plus the work rate due to pressure forces is 
written as 


a 2 (3 73 ~2 h 2 


~ CH ~37 CT $ 7 2 —3 7 ~7 ~72 

qe = —u h -\ - —b%u + abuh + auh 


0?P7.3~( = 


h u[ u~j + auuxx - -cm- — 


a 3 fir - 


I) 


bxUxU h 


272 


abxx'U 2 + h (auux + u t ~)). (52) 
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With these definitions, the energy balance is 




( 53 ) 


Using the scaling E = pc^b^E and qE — PCq&o9e, the dimensional form of energy 
density per unit span in the transverse direction is given as the sum of the kinetic 
and the potential energy by 

E = ^u 2 (l + b 2 x )h - ^ uu x b x h 2 + ^ 2 h 3 + ^-h 2 + pgbh , (54) 

^--V----' ^-v-' 

Ek E p 


and the dimensional form of energy flux plus work rate due to the pressure force is 
given by 


Qe = pgu{h + bh) 


i 3 h{l + b 2 x ) - ^h 3 u(^- 


u[u xt + uu xx - -ul^ 


P- 2- 7 j2 I 

2 'U' 'U'xOxh' H - 


i ah i 


b'r.'r'U 


b x (i 


U t 


(55) 


For a horizontal bed, it is more convenient to normalize the potential energy of a 
fluid particle to be zero at the bottom. If this is done, then the dimensional forms 
of energy density and energy flux plus work rate due to pressure forces are given 

by 

E = f A 2 + + Zh>«i , (56) 

and 


- 2 P -3 P 3-(- _ 3 _2\ 

Qe pguh —u h — h uyuxt H - 'U'Uxx ~^2^ x ) ’ 


(57) 


respectively. Note that as (3 —> 0 in the equation (l49l) . the energy balance reduces 
to the shallow-water energy conservation with 


£ s “ = f (b 2 +2b 0 r 1 + V 2 ) + ^hu 2 


(58) 


and 


q S E = pgh 2 u + ^h.u 3 . (59) 

In addition, in the case a ~ /?, the energy balance reduces correctly to the case 
of the classical Boussinesq system, with E = ^(&q + 2bog + rj 2 ) + %b 0 u 2 and 
Qe = pg(bl + 2gb 0 )u. 

It is worth noting that the conservation of the asymptotic approximation to the 
total energy with nontrivial bathymetry in the fully nonlinear regime is satisfied by 
the solutions of the SGN equations exactly. This can be seen by performing lengthy 
computations using formal integrations by parts, or by recognizing that the total 
energy £ = J E dx differs from the Hamiltonian dll) by a term | f b 2 dx which is 
constant. 
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4. Applications 

4.1. Evolution of undular bores. In free surface flow, the transition between 
two states of different flow depth is called a hydraulic jump if the transition region 
is stationary, and a bore if it is moving. Bores are routinely generated by tidal 
forces in several rivers around the world, and may also be generated in wavetank 
experiments Kj. ;44j. 

The experimental studies of [43] show that when the ratio between the differ¬ 
ence in flow depths to the undisturbed water depth is smaller than 0.28, then the 
bore will feature oscillations in the downstream part. If this ratio is greater than 
approximately 0.75, then a so-called turbulent bore ensues. If the ratio is between 
0.28 and 0.75, the bore will be partially turbulent, but will also feature some os¬ 
cillations. The bore strength can also be expressed in terms of the Froude number 
Fr = \J [(2/ii //iq + l) 2 — l]/8, and more recent studies, such as [44] have found that 
when Fr > 1.4 approximately, the bore consists of a steep front, while undulations 
are growing at the bore front only in the near-critical state Fr ss 1. The different 
shapes and a transition from the subcritical to the supercritical regime is described 
in [45], and in [46] an empirical critical value Fr cr it = 1-3 is suggested in order to 
determine the breaking of an undular bore. However the exact characterization of 
the transition between these states still remains unclear. 

Some of the divergence in the results on the critical bore strength might be 
explained by the observation that one single nondimensional number may not be 
sufficient to classify all bores. For example, in [47], a hyperbolic shear-flow model 
is suggested which allows the classification of bores with an additional parameter 
depending on the strength of the developing shear flow near the bore front. 

The connection between the initial bore strength and the ensuing highest un¬ 
dulation is fairly well understood. Using Whitham modulation theory [48] . it can 
be shown that if viscosity is neglected, the amplitude of the leading wave behind 
the bore front is exactly twice the initial ratio of flow depths [41 ] 149] . This result 
agrees well with experimental findings. For example, the amplitude of the leading 
wave found experimentally in [43] was 2.06 times the initial amplitude ratio. 

In this section we present a numerical study of the energy balance of undular 
bores for the SGN equations. It is well known that a sharp transition in both flow 
depth and flow velocity which conserves both mass and momentum necessitates 
a loss of energy across the front. The standard argument essentially relies on an 
inviscid shallow-water theory and the examination of exact weak solutions of the 
shallow-water equations [29]. 

Given these assumptions, it is natural to explain the energy loss across the bore 
front by pointing to the physical effects neglected in the shallow-water theory, such 
as viscosity, frequency dispersion, and turbulent flow. Indeed, in strong bores, 
turbulent dissipation accounts for the lion’s share of energy dissipation, and a long¬ 
wave model can only give a first approximation of the dynamics. Most of the work 
investigating the energy loss has focused on weak undular bores, where long-wave 
models can be expected to yield an accurate description of the flow. The loss of 
energy in weak bores has been explained by the creation of oscillations in the free 
surface behind the front, but it was noted in [30] that an additional dissipation 
mechanism is needed. In [31] . the bottom boundary layer was invoked to explain 
this required additional energy loss, but it was noted in 1321 33] that invoking 
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frictional effects to explain the energy loss experienced by a conservative system 
was not consistent. 

However, as already mentioned, there was a slight technical problem in the anal¬ 
ysis of [33j . since the energy functional 


£b 


ous 



a 2 w 2 h + (wwxx + + h 2 


dx , 


(60) 


used in that work could not be obtained in the framework of the asymptotically 
correct mechanical balance laws derived in [55]. Indeed, the expressions for the 
energy and energy flux associated to the Boussinesq system which were derived in 
[26] are 


E = 


arj - 


ct 2 ~2 


a 2 _2 a 12 l G;2 ~2 

T w = 2 h + T” 


and the non-dimensional energy flux (corrected for the work rate due to pressure 
forces) as 


_ 2 -- a P 

qE = olw + 2a wrj H—— 


1 


afd 


d - - ) Wxx = cm + 2a wrj + —wss , 


where w is the nondimensional horizontal velocity component at height 6 = 6 0 -^/2/3 
in the water column. It is apparent that as /3 —> 0, these expressions do not reduce 
correctly to the corresponding expressions of the shallow-water theory. However, 
since the expressions (1551) and do reduce to the correct shallow-water equiv¬ 
alents, the analysis of the energy loss in the undular bore can be made precise in 
the context of the SGN system. 

The numerical method that was used to perform the numerical simulations in 
this paper is detailed in El It is also noted that for simplicity’s sake we consider 
the water density p — 1 kg/m 3 . The numerical experiments require initial data. 
An initial surface condition that triggers the generation of undular bores is 

h(x,0) = ho + -(/li — ho)tanh(ra) , 

where k is the parameter that determines the steepness of the undular bore. Here 
we take n = 1/2. In order to generate a simple undular bore, i.e. a wave that 
propagates mainly in one direction, we consider an initial flow given by the following 
velocity profile: 


Sh ( \ 3 ^ 2 

u(x, 0) = — ( (2 /iq + 3(5h)ho + (dh) 2 ) ) x (1 — tanh(rea:)) , 

hi \2ho / 

where Sh = hi — ho- One may envision other numerical methods to create an 
undular bore, such as the addition of a line source in the upstream part, such as 
used in [50]. Nevertheless, the initial conditions described above were sufficient for 
our purposes. 

First we present the computation of the energy budget in an undular bore for 
various bore strengths. We consider the control volume [xi,X 2 ], where x\ is far to 
the left of the bore front, and X 2 is far to to the right. In Table [T] the bore strength 
is shown in the first column, and the corresponding Froude number is shown in 
the second column. Taking ho = 1 and h\ between 1.1 and 1.7 we monitor the 
energy flux and work rate due to the pressure force, given by qe{x\) — qE{x 2 ) as 
defined in <5ZD, and these values are shown in the third column of the table. We 
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Table 1. Energy conservation 


hi/ho 

Fr 

qE{xi) - qs{x 2 ) 

d£/dt 

1.1 

1.07 

3.6481059 

3.6481059 

1.2 

1.15 

8.6017456 

8.6017456 

1.3 

1.22 

15.100378 

15.100378 

1.4 

1.30 

23.394470 

23.394470 

1.5 

1.37 

33.746103 

33.746103 

1.6 

1.44 

46.429376 

46.429376 

1.7 

1.51 

61.730669 

61.730669 

Table 2. 

Momentum conservation 

hi/ho 

Fr 

qi(x 1 ) - qi{x 2 ) 

dX/dt 

1.1 

1.07 

1.1330550 

1.1330549 

1.2 

1.15 

2.5898340 

2.5898399 

1.3 

1.22 

4.3997850 

4.3997849 

1.4 

1.30 

6.5923199 

6.5923198 

1.5 

1.37 

9.1968749 

9.1968748 

1.6 

1.44 

12.242880 

12.242879 

1.7 

1.51 

15.759764 

15.759764 


also monitor the gain in energy in the control interval as given by £(t) = f '^ E dx. 
These values are shown in the fourth column. The particular figures shown in the 
table are for T = 30, but the values are nearly constant over time. It is apparent 
from the table that energy conservation holds to at least eight digits, even for large 
bore strengths. These numbers confirm our previous finding that the energy is 
exactly conserved in the SGN model, and also validates the implementation of the 
numerical method. In addition, these results confirm our claim that no dissipation 
mechanism is necessary to explain the energy loss in an undular bore. 

As noted in the previous section, the expression (1571) for the energy flux and work 
rate due to pressure forces reduces to the corresponding formula for the shallow- 
water theory in the case of very long waves. Since x\ and X2 are relatively far from 
the bore front, shallow-water theory should be valid at these points. Therefore, the 
usual formula for the energy loss in an undular bore in the shallow-water theory is 
valid: 

rfTpsw n I Z T” 

+q S E(x2,t) - q S E(x!,t) = -|(/h - h 0 ) 3 )J Iff 3 + 7 ^-J • (61) 

Since there is no energy loss in a dispersive system, one may conclude that the excess 
energy is fed into oscillations of the free surface, and the formula (1611) furnishes an 
estimate of the amount of energy which is residing in the oscillatory motion. 

A similar study can be performed on the momentum balance. Momentum gain in 
the control interval is given by the momentum flux through the lateral boundaries 
and the pressure force as qi{x\) — qi(x 2 ), with qi given in (l43l) up to T = 30. Table 
[2] presents the momentum rates. As in the case of the energy, the corresponding 
values agree to about eight digits. In Figure [5] we present the normalized values 
0 ) of the momentum and £{£)/£(0 ) of the total energy for the values of the 
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Figure 2 . The momentum and the energy of the undular bore for 
Fr = 1.07, 1.30 and 1.51. 


Froude number Fr = 1.07, 1.30 and 1.51. The slopes of the lines can be found in 
Tables |T| and [2] 

Figure [3] shows the profiles of the undular bores generated when hi/ho = 1.3, 1.5 
and 1.7. From these figures we observe that as the Froude number Fr increases, 
the peak amplitude of the leading wave becomes larger, and the shape of the wave 
envelope is changing. For example the shape of the wave envelope of Figure [3] (a) 
can be described by a quadratic function (wineglass shape) while the shape of the 
wave envelope of Figure [3] (c) can be described by a square-root function (martini 
glass shape). For the various shapes of the undular bores we refer to [511 . 


4.2. Shoaling of solitary waves. In this section we study the conservation of 
energy in the case of a nonuniform bathymetry. Specifically, we consider the ex¬ 
periments proposed in [551 [55] related to the shoaling of solitary waves on a beach 
of slope 1 : 35. The shoaling of solitary waves has been studied theoretically and 
experimentally in many works, such as in [52j[53] HU[55]. Next, we study the shoal¬ 
ing of solitary waves with normalized amplitude A = 0.1, 0.15, 0.2 and 0.25 in the 
domain [—100,34]. In the numerical experiments we take Aa; = 0.05 while we trans¬ 
late the solitary waves such that the peak amplitude is achieved at x = —20.1171 
while the bottom is described by the function 


f —1, x < 0 

( —1 + a;/35, x > 0 ’ 
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Figure 3. Undular bores profiles for various Fr values. 


but modified appropriately around x = 0 so as to be smooth enough and to satisfy 
the regularity requirements of the model. 

A comparison between the experimental results on shoaling waves from [52J, and 
the shoaling solitary waves computed with numerically approximation of © and 
(El is presented in Figured Overall, we observe a very good agreement between 
the numerical results and the experimental data. 

Tabledpresents the conserved values of the total energy £ and of the Hamiltonian 
H for t £ [0, 45] for the computations shown in Figured We observe that the energy 
is conserved with more than ten decimal digits. Due to the small values of Ax and 
At no energy dissipation can be observed verifying the efficacy of the numerical 
method. 


Table 3. Conserved values of energy (in Joules) and Hamiltonian 
for shoaling of solitary waves on a plane beach of slope f : 35. 


A 

£ 

n 

0.10 

62.4704607870 

0.05202930490 

0.15 

62.7102258381 

0.09856973753 

0.20 

62.9401348199 

0.15627417412 

0.25 

63.1680884219 

0.22460417742 























18 


HENRIK KALISCH, ZAHRA KHORSAND, AND DIMITRIOS MITSOTAKIS 



Figure 4. Comparison of the numerical solution and the experi¬ 
mental data on wave gauges of [53]. — Numerical solution;-: 

Experimental data. 


Although the total energy is conserved the kinetic and the potential energy are 
not constant with time. Figure [5] presents the normalized kinetic energy £k(t)/£k{ 0) 
and normalized potential energy £ p (t)/£ p ( 0) evaluated in the spatial interval [—100,34]. 
As can be seen in Figure [5] the kinetic energy is decreasing at the early stages of 
shoaling due to the slight decrease in the wave speed while the potential energy 
is initially increasing due to the increase of the wave height. At later stages of 
the shoaling, the kinetic energy increases again, due to the increase in particle ve¬ 
locities, and the potential energy decreases again, due to the rising bottom, and 
narrowing wave peak. Nevertheless, the total energy is constant over time. 

5. Summary and Conclusions 

We have detailed the derivation of mechanical balance laws for the SGN equa¬ 
tions in the case of a horizontal bed and also in the case of varying bathymetry. The 
mechanical balance laws derived here, including the mass, momentum and energy 
balance laws, are valid to the same asymptotic order as the SGN system, providing 
a firm link between conservation laws associated to the governing SGN equations, 
and the above mechanical quantities. Finally, applications to the energy budget 
of undular bores and the development of potential and kinetic energy in shoaling 
solitary waves have been presented. In particular, it has been shown that the en¬ 
ergy loss in undular bores is fully compensated for by the development of surface 
oscillations, since the energy in the SGN with a flat bottom is exactly conserved. 
Indeed, exact conservation of energy to near machine precision was observed in our 
numerical computations, and this gave an additional check on the implementation 
of the numerical algorithm. 
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Figure 5. Normalized Kinetic and Potential energy for shoaling 
of solitary waves on a plane beach of slope 1 : 35. 
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Appendix A. The numerical method 


In this Appendix we consider the initial-boundary value problem (IBVP) com¬ 
prised of system (I28al) (I28bl) subject to reflective boundary conditions. Rewriting 
the system in terms of (h, it), and dropping the bar over the symbol of the hori¬ 
zontal velocity, yields the IBVP 


h t + ( hu) x = 0 , 

hu t + huu x + g(h + b) x + 


[h 2 a'P+lQ)] x + b x (lV+Q)= 0 , 

u(a, t ) = u(b, t) = 0 , 
h[x, 0 ) = h 0 (x) , 
u(x,0) = u 0 (x) , 


(62) 


where V = h [u 2 — u xt — uu xx \ , Q = b x (u t + uu x ) + b xx u 2 , x € [a, i] Cl and 
t € [0, T], Considering a spatial grid Xi = a + i Ax, for i = 0,1, • • • , N, where Ax 
is the spatial mesh-length, such that Ax = (6 — a)/N, JVeN. We define the space 
of cubic splines 


S = {</> € C 2 [a,b\ <A|[ Xi , Xi+l] € P 3 , 0 < i < N - l| 
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where P fc is the space of polynomials of degree k. We also consider the space 
= S n {</> £ C([a, b ]) (j)(a) = = o} . 


The basis functions of the space S and S 0 consist of the usual B-splines described 
in [55] . 


N 

E 0 [H\ 

rate for E 0 [H] 

E 0 [U\ 

rate for E 0 [U] 

300 

0.1211 x 10 -s 

- 

0.6127 x 10" 11 

- 

320 

0.9674 x 10~ 9 

3.4793 

0.4733 x 10" 11 

3.9983 

340 

0.7836 x 10" 9 

3.4772 

0.3714 x 10 -11 

3.9999 

360 

0.6422 x 10~ 9 

3.4797 

0.2955 x 10 -11 

3.9977 

380 

0.5322 x 10 -9 

3.4754 

0.2382 x 10" 11 

3.9885 

400 

0.4452 x 10 -9 

3.4793 

0.1939 x 10~ n 

4.0099 


Table 4. Spatial errors and rates of convergence in the L 2 norm. 


The semi-discrete scheme is reduced in finding h £ S and u £ So such that 
(h t ,<p) + {ihu) x ,(^ = 0 , 

B(u t ,ip] h) + (Jiuu x + g(h + b) x , ip'j + ( 63 ) 

+ (h 2 (lV + ±Q),^) + (b x (\V + Q),i) =0 , 

for (f> £ S, and ip £ So, and V = h [u 2 — uu xx ] and Q = b x uu x + b xx u 2 . B is defined 
as the bilinear form that for fixed h is given by 


= h 


1 h x b x T ^ hb xx T b x 


+ 5 ( h 3 ip x ,Xx ) for ip,x £ S 0 ■ 


The system of equations (1631) is accompanied by the initial conditions 
h{x , 0) = V{h 0 (x)} , u(x, 0) = Vo{uo{x)} , 


(64) 

(65) 


where V and Vo are the L 2 -projections onto S and So respectively, satisfying 
(Vv,(p) = {v,(f>) for all 0 £ S and (Vov,ip) = (v,ip) for all ip £ So- Upon choosing 
basis functions <pj and ipj for the spaces S and So, (l63l) is reduced to a system of or¬ 
dinary differential equations (ODEs). For the integration in time of this system we 
employ the Dormand-Prince adaptive time-stepping methods, mm- One may 
apply the same numerical method to solve the IB VP with non-homogeneous Dirich- 
let boundary conditions. For example if u(a , t) = u a then the change of variables 
u(x,t) = w(x,t) + uq(x) reduces the non-homogeneous system to a homogeneous 
IB VP system for the variable w. In all the numerical experiments we took Ax = 0 . 1 , 
while the tolerance for the relative error of the adaptive Runge-Kutta scheme was 
taken 5 • 10 -14 . For the computations of the integrals the Gauss-Legendre quadra¬ 
ture rule with 8 nodes was employed. 

The convergence properties of the standard Galerkin method for the SGN system 
are very similar with those of the classical Boussinesq system studied in detail in 
|59l 160 ] . In order to compute the convergence rates in various norms we consider 
the nonhomogeneous SGN system with flat bottom admitting the exact solution 
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h(x,t) = 1 + e 2t (cos( 7 ra;) + x + 2) and u(x,t) = e tx x sin(7rx) for 0 < x < 1, and 
for t £ (0,T] with T = 1. We compute the normalized errors 


p rpi II F{x, T', Ax) Fexa,ct( x i TIU 

Sl J ITexactMIU 


( 66 ) 


where F = F(-;Ax) is the computed solution, i.e., either H « h{x,T) or U ~ 
u(x,T), inexact is the corresponding exact solution and s = 0,1, 2, oo correspond to 
the L 2 , H 1 , H 2 and L°° norms, respectively. The analogous rates of convergence 
are defined as 


= ln(£' a [F(-; Axk-i)]/E s [F(-\ Axfc)]) 

ln(Aa; fc _i/Axfc) 

where Axk is the grid size listed in row k in table [J] To ensure that the errors 
incurred by the temporal integration do not affect the rates of convergence we use 
At < Ax while we take Ax = 1/N. 

Table [4] presents the spatial convergence rates in the L 2 norm. We observe that 
the convergence is optimal for the u variable but suboptimal for the h variable. 
Specifically, it appears that \\h— h|| ~ Ax 3 - 5 , while ||u — u|| ~ Ax 4 . More precisely, 
as in the case of the classical Boussinesq system J>5], and because the rate of 
convergence in h appears to be less that 3.5 yields that the error should be of 
0(Aa; 3 ' 5 ^/ln(l/A.T)). Similar results obtained for the convergence in the H 1 , H 2 
and L°° norms. Specifically it was observed numerically that \\h — h\\ s ~ Acc 3 ' 5_s , 
||w — w||i ~ Ax 4_s , for s = 0,1, 2 and \\h — h||oo ~ Aa; 3 , while ||w — uHoo ~ Ax 4 
approximately. 


rate for E S [F] 
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